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Genetic robustness, the preservation of an optimal phenotype in the face of mutations, is critical 
to the understanding of evolution as phenotypically expressed genetic variation is the fuel of natural 
selection. The origin of genetic robustness, whether it evolves directly by natural selection or it is 
a correlated byproduct of other phenotypic traits, is, however, unresolved. Examining microRNA 
(miRNA) genes of several eukaryotic species, Borenstein and Ruppin (Borenstein et al. 2006, PNAS 
103: 6593), showed that the structure of miRNA precursor stem-loops exhibits significantly increased 
mutational robustness in comparison with a sample of random RNA sequences with the same stem- 
loop structure. The observed robustness was found to be uncorrelated with traditional measures 
of environmental robustness - implying that miRNA sequences show evidence of the direct evo- 
lution of genetic robustness. These findings are surprising as theoretical results indicate that the 
direct evolution of robustness requires high mutation rates and/or large effective population sizes 
only found among RNA viruses, not multicellular eukaryotes. We demonstrate that the sampling 
method used by Borenstein and Ruppin introduced significant bias that lead to an overestimation of 
robustness. Introducing a novel measure of environmental robustness based on the equilibrium ther- 
modynamic ensemble of secondary structures of the miRNA precursor sequences we demonstrate 
that the biophysics of RNA folding, induces a high level of correlation between genetic (mutational) 
and environmental (thermodynamic) robustness, as expected from the theory of plastogenetic con- 
gruence introduced by Ancel and Fontana (Ancel et al. 2000, J. Exp. Zool. 288: 242). In light of 
theoretical considerations we believe that this correlation strongly suggests that genetic robustness 
observed in miRNA sequences is the byproduct of selection for environmental robustness. 



Introduction 



The magnitude of genetic effects on phenotype depends 
strongly on genetic background, the effects of the same 
mutation can be larger in one genetic background and 
smaller in another. The idea that wild-type genotypes 
are mutationally robust, i.e. show invariance in the face 
of mutations (more gene rally heritable per turbations), 
goes back to Waddington Waddington 1957j , who origi- 
nally introduced the concept as canalization. While ge- 
netic robustness has been found across different levels of 
organization from individual genes, through simple ge- 
netic circuits to entire organisms (approximately 80% 
of yeast single knockouts have n o obvious effect in rich 
medium Hillenmever et al. 2008| |). the origin of the ob- 
served robustness has remained a source of contention. 
The three main hypotheses regarding the potential origin 
of genetic robustness predate the concept itself, and fall 
along the lines of the famous debate between members of 
the modern synthesis (in particular Wright, Haldane and 
Fisher) surrounding the origin of dominance (dominance 
can be understood as a simple case of robustness, a dom- 
inant phenotype being more robust against mutations) 
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de Visser et al. 2003L iMavo and Burger 19971 ]: (i) the 
most straightforward explanation, favoured by Wright, 
was tha t robustness evolves directly, through natural se- 
lection Fisher 19281 ]: (ii) an alternative congruent hy- 
potheses, put forward in the context of dominance by 
Haldane, proposes that the evolution of genetic robust- 
ness is a correlated byproduct of selection for environ- 
mental robustness, i.e. invariance in the face of nonher- 
itable perturbations, e.g. temperature, salinity or inter- 
nal factors such as fluctuations in the concentrat ion of 
gene products during development [Haldane 1930j ; (iii) 
while a third view holds that genetic robustness is in- 
trinsic, arising simply because the buffering of a charac- 
ter with respect to mutations is the necessary or likely 
consequence of cha racter adapta tion, in the context of 
domin ance Wright Wright 19341 and later Kacser and 
Burns [Kacser and Burns 1981 1 argued that it arises as 
an inevitable, passive consequence of enzyme biochem- 
istry and selection for increased metabolic flux. 

Recently robustness has been a subject of re- 
newed interest. Several theoretical and simula- 
tion studies have addressed robustness in a wide 
range of contexts rangin g from gene redundancy 
Krakauer and Plotkin 20021 to model regulatory net 
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works Siegal and Bergman 20021 Azevado et al. 2007 . 
Ciliberti et al. 2007al. ICiliberti et al. 2007b . 

Crombach and Hogeweg 20081 ]. By in the first 
case building on evidence of excess mutational 
robustness present in RNA secondary structure 



Wagner and Stadler 1999| and in the second case 
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on the expectation that high mutation rates present 
amo ng RNA viruses should favour mutational robust- 
ness (Wilke and Adami 2003ll two pione ering studies, by 
Montville et al. Monetville et al. 2005 1 an d Borenstein 
and Ruppin Borenstein and Ruppin 20061 ] have man- 
aged to step beyond computer simulations and through 
using, respectively, i n vitro evolution experiments 
Monetville et al. 2005| and microRNA sequences from 



diverse taxa found evidence to support the hypotheses 
that genetic robustness can evolve directly. Further 
work on in vitro evolution experiments has provided 
additional evidence showing that if a population is 
highly polymorphic robustness can evolve directly 
jSanjuan et al. 2007L [Bloom et al. 200T\ . 

The theoretical underpinnings of these studies is 
provided by the results of van Nimwegen et al. 
van Nimwegen et al. 19991 ] . who through solving the 
quasispecies equations describing the evolution of a 
population on a network of phenotypically neutral se- 
quences, were able to demonstrate, that provided a 
sufficiently polymorphic population, mutational robust- 
ness can evolve directly. The necessary mutation rates 
and/or population sizes were found to be very large in 
simulation studies using R NA secondary structure as 
a genotype-phenoty p e map |van Nimwegen et al. 19991 
iForester et al. 20061 . ISzdllosi and Derenvi 20081 ] . direct 
evolution of increased neutrality requiring the product 
of the effective population size N e and the mutation rate 
per nucleotide u to be well in excess of one. Such high 
mutation rates can only readily be found among RNA 
viruses, are extraordinary even among unicellular organ- 
isms (Prochlorococcus 2N e u f=a 2., E. coli 2N e u ~ 0.2, 
S. cerevisiae AN e u f=a 0.09) and completely unheard 
of among multicellular eukaryotes possessing RNA si- 
lencing mechanisms and microRNA genes (A thaliana 
AN e u « 0.012, D. melanogaster 4N e u w 0.015, C. ele- 
gans 4N e u w 0.013, C. intestinalis AN e u w 0.012, M. 
musculus AN e u « 0.001, H. sapiens AN e u w 0.001) 



Lynch and Conerv 20031 ] . 



In their study Borenstein and Ruppin examined mi- 
croRNA (miRNA) precursor sequences from several eu- 
karyotic species. miRNA arc small endogenous non- 
coding RNAs that regulate the expression of protein 
coding ge nes through the RNA int er ference (RNAi) 
pathway lLagos-Quintana et al. jOOll E au et al. 20011 



iLee and Ambros 200l[ iBartel 20041 ] . Functionally rele 
vant short (« 22nt) mature miRNA sequences are ex- 
cised from longer precursor sequences that fold into a 
stem-loop hairpin structure. The hairpin like secondary 
structure of precursor s tem-loops pl ays a crucial role in 
the maturation process Bartel 2004j and is under evolu- 
tionary constraint to conserve its structure. Borenstein 
and Ruppin used the novel and ingenuous method of 
generating for each miRNA sequence a random sample 
of sequences with identical minimum free-energy (MFE) 
structure to uncover traces of adaptation. To compare 
the mutational robustness of miRNA precursor sequences 
to random sample sequences with identical MFE struc- 



ture they compared the single mutant neighborhood of 
a given miRNA precursor sequence to the single mutant 
neighborhood of the sample sequences. Calculating the 
average distance of the MFE structure of each single mu- 
tant sequence to the MFE structure of the original se- 
quence for both stem-loop and sample sequences (details 
on secondary structure calculations are presented be- 
low) they demonstrated that miRNA precursor sequences 
have single mutant neighborhoods with sequences that 
fold into more similar MFE structures compared to se- 
quences in the single mutant neighborhoods of sample 
sequences with identical MFE structure. While a similar 
comparison of the folding minimum free-energy showed a 
comparable, but lower bias, the finding that the two were 
only weakly correlated allowed the authors to conclude 
that the observed bias is a result of direct selection for 
mutational ro bustness. Their results were reexamined 
by Shu et al. Shu et al. 2007| who argued that muta- 
tional robustness among miRNA precursors may be the 
correlated byproduct of selection for environmental ro- 
bustness, but found only a moderately higher correlation 
using a different measure of mutational robustness. 

In light of the consistently low value of uN e among 
multicellular eukaryotes the results of Borenstein and 
Ruppin are highly surprising. There is no known mecha- 
nism which can explain the direct evolution of robust- 
ness that they observe. According to the classic re- 
sults of Kimura and Maruyama the average fitness of 
an asexually reproducing population (in the limit of very 
large populations) depends only on the mutation rate 
and is independent of the details of the fitness land- 
scape [Kimura 19661 ]. This result, however, only holds 
under the assumption that the fittest genotype does not 
have any neutral sites. While, the extension of these 
results to more general fitness landscapes by van Nimwe- 
gen et al. demonstrates that the presence of neutral 
genotypes can lead to selective pressure to evolve mu- 
tational robustness simulation studies using genotype- 
phenotype maps induced by RNA secondary structure 
have demonstrated that uN e > 1 is a necessary condi- 
tion Forester et al. 20061 even n th e presence of rccom- 
bination [Szollosi and Derenvi 2008| . The case for the 
direct evolution of genetic robustness rests on the paral- 
lel findings that a stronger bias for mutational robustness 
is present in miRNA precursor sequences than for envi- 
ronmental robustness and that the two are only weakly 
correlated. Introducing a new measure of environmen- 
tal robustness in this paper we endeavor to demonstrate 
that, indeed as previously also suggested by Shu et al. 
Shu et al. 2007| |. the exact opposite is true: the bias 



for environmental robustness is stronger and it is highly 
correlated with mutational robustness. The correlated 
evolution of environmental and mutational robustness in 
RNA sequences under selection to retain secondary struc- 
ture is expected as a corollary of a general casual link 
between environmental robustness and genetic robust- 
ness in RNA sequences p roposed by Ancel and Fontana 
' Ancel and Fontana 2000 ]. They argue that plastogenetic 
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enna RNA package [Hofbacker et al. 1994 1 to produce a 
sequence with MFE structure identical to that of the na- 
tive sequence that is stored (ii) and subsequently random- 
izing this sequence by attempting 4L random nucleotide 
substitutions in a miRNA precursor sequence of length L, 
accepting a substitution if the resulting sequence's MFE 
structure remains unchanged. For each miRNA precursor 
sequence on average > 800 sample sequences with identi- 
cal MFE structure was generated. Supplemental material 
accompanying our paper contains the robustness values 
for all 4361 genes associated with 3641 unique sequences 
we considered. 



FIG. 1: a Generating a random sample of sequences 
with a desired MFE structure by stochastic minimization 
of the free-energy of the desired fold (the method em- 
ployed t he RNAinverse program use d by Borenstein and 
Ruppin iBorenstein and Ruppin 20061 ] as well as Shu et al. 
[Shu et al. 2007t | ) results in a biased sample in which se- 
quences with lower than average neutrality (higher than av- 
erage number single mutant neighbors) are overrepresented. 
This can be avoided if after finding a sequence with the de- 
sired MFE structure a random walk is performed among se- 
quences with the desired MFE structure. This random walk 
on the neutral network associated with the MFE structure 
mimics the sequnce drift of a sequence evolving under the 
constraint to fold in to the desired MFE structure, b Rank 
score distributions for two measures of mutational robust- 
ness {r} s and r\ n see text). Comparing the distributions de- 
rived from sampling using only stochastic optimization (top, 

f bia S ed = q ^biased = q f biascd = q 3^ ^biased = Q gg) 

to that derived from sampling with subsequent randomiza- 
tion (bottom, f a = 0.29, R B = 0.78, f n = 0.44, R n = 0.59) 
shows that increased neutrality is predominately an artifact 
of biased sampling, while the lower than average distance of 
MFE structures in the mutational neighborhood to the wild- 
type MFE structure becomes somewhat less pronounced, but 
is still significant. 



congruence, i.e. the correlation between the set of struc- 
tures thermally accessible to a sequence, its plastic reper- 
toire, and the MFE structures of its genetic neighbour- 
hood will lead to the emergence of mutational robustness 
in the presence of selection for some predefined structure. 



Materials and Methods 
microRNA sequences and sampling 

miRNA precursor s equences were downloaded from 
miRBase version 9.0 [Griffiths- Jones et al. 200"f| . All 
4361 miRNA genes were used, yielding 3641 unique 
miRNA precursor sequences. For each miRNA precursor 
sequence we produced a sample of random sequences by 
(i) using the stochastic optimization routine from the Vi- 



Measuring thermodynamic robustness 

In order to calculate the thermodynamic robustness 
measure r\ t , defined below, we sampled the equilibrium 
thermodynamic ensemble of stem-loop and sample se- 
quences using the stochastic backtracking routine from 
the Vienna RNA package producing 10 6 suboptimal 
structures per sequence, using the default temperature 
of 310 K. The average distance from the MFE struc- 
ture in the thermodynamic ensemble can be calculated 
exactly with the help of base-pairing probabilities, which 
are available as a byproduct of partition function fold- 
ing in the Vienna package, and were used to validate the 
sampling. 



Statistics 

Given a rank score r and sample size N a good estimate 
for the probability of observing an equal or lower rank 
scor e by chance is given by (rN) /(N + 1) w r. Follow- 
ing [Borenstein and Ruppin 20061 ] rank scores of r < 0.05 
are considered significantly robust. To determine if the 
robustness of miRNA precursor sequences according to 
some measure rj has the same distribution as the robust- 
ness of s ample sequences r\' for a gro up of sequences, 
following [Borenstein and Ruppin 20061 ] we test against 
the null hypothesis that they are drawn from identical 
distributions using the n onparametric Wilcoxon signed 
rank test. In contrast to Borenstein and Ruppin 20061 ] . 
however, we do not consider as paired values r\ and the 
average of rf over all N sample sequences, (77'), as we 
found this to result in spuriously low p- values, but instead 
calculate the p- values for a given group of sequences by 
averaging over 1000 different sets of {77, 7/} pairs where 
in each set the rf values belong to a random sample se- 
quence. As a complementary approach we also tested 
the hypothesis that the distribution of rank scores of a 
group of sequences for a given robustness measure is uni- 
form - as we would expect if miRNA precursor sequences 
were randomly sampled from the set of sequences with 
identical MFE structure - using a standard Kolmogorov- 
Smirnov goodness of fit test. We found the two signifi- 
cance analyses to be in good agreement indicating highly 
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significant bias for higher values of r]t{dth.) an d Vs, but 
mostly no or only nonsignificant bias for higher r) n . The 
supplemental information accompanying our paper con- 
tains species level statistics and significance analyses. 



Results 

We assessed the environmental and mutational robust- 
ness of 3641 unique miRNA precursor sequences and for 
each sequence compared them to a random sample of 
sequences with the same MFE structure. The idea of 
looking for signs of adaptation for increased robustness 
among miRNA precursor sequences by comparing the ro- 
bustness of naturally occurring sequences to that of ran- 
dom sequences with the same secondary structure is con- 
ceptually similar to the approach used to support the 
argument that t he genetic code has evolved to minimize 
mutational load |Pi Giulio 19871 lHaigh and Hurst 199lT 
ISzathmarv and Smith 1997 1. In the case of the genetic 
code the authors took the common genetic code, and, for 
each codon, calculated the change in polarity of the en- 
coded amino acid caused by replacing each of the three 
nucleotides, one after the other. In order to determine 
whether the genetic code is adapted to minimize mu- 
tational load they proceeded by comparing the mean 
squared change caused by the replacement of a single nu- 
cleotide in the common genetic code to 10000 randomly 
generated codes with the same redundancies. They found 
that only two of the random codes were more conser- 
vative than the common code with respect to polarity 
distances between neighbouring amino acids. 

We undertook a similar program in the case of miRNA 
precursor sequences. Each miRNA gene encodes a short 

22 nucleotide sequence that is partially complemen- 
tary to the mRNA of proteins regulated by the partic- 
ular miRNA gene. For the proper short sequence to be 
excised by the protein Dicer, and hence for the miRNA 
gene to be functional, a larger part of the miRNA se- 
quence, called the miRNA precursor sequence, must fold 
into the proper secondary structure. In order to deter- 
mine whether a miRNA precursor sequence is adapted to 
minimize the effects of mutational and/or environmental 
perturbations, i.e. to maximize mutational and/or envi- 
ronmental robustness, we compared the mutational and 
environmental robustness of each miRNA precursor se- 
quence (robustness measures we used are defined below) 
to the mutational and environmental robustness of a ran- 
dom sample of sequences with identical structural phc- 
notype (i.e. identical MFE structure). 

To generate a random sample of sequences with given 
MFE structure we first used, starting from a random 
sequence, stochastic minimization of the free-energy of 
the target structure to find a sequence with the de- 
sired MFE structure. This method by itself, however, 
yields a biased sample of sequences (see Fig. QJi and 
b) and must be supplemented by an additional ran- 
domization step (see Materials and Methods). To mea- 
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FIG. 2: In order to examin the robustness of miRNA precur- 
sor sequences to thermal fluctuations we sampled the equilib- 
rium thermodynamic ensemble of structures. Sampling 10 6 
structures for each miRNA precursor sequence and each mem- 
ber of the random sequence sample, we binned structures ac- 
cording to their distance from the MFE structure, a For, e.g. 
the Monodelphis domestica miRNA precursor sequence mdo- 
mir-1 examining the distribution of structures as a function 
of the base-pair distance shows that the averaged random se- 
quence sample distribution (white bars) has a much larger 
fraction of structures that are drastically different from the 
MFE structure, compared to the distribution of structures for 
the original miRNA precursor sequence (black bars), b Ex- 
amining the averaged distribution of stem-loop (black bars) 
and random corresponding random sequence sample distri- 
butions (white bars) shows that there is a general tendency 
among miRNA precursor sequences for increased thermody- 
namic robustness, i.e. of avoiding structures that are highly 
dissimilar to the MFE structure. A strikingly similar effect 
can be observed if we examine the distribution of structures in 
the mutational neighborhood. Analogous to a, in c we binned, 
according to their distance from the MFE structure of the wild 
type, the MFE structures of all (3L) single point mutants for 
the Monodelphis domestica miRNA precursor sequence mdo- 
mir-1 (black bars) as well as the MFE structures for each se- 
quence in the single mutant neighborhood of sample sequences 
(white bars). The distribution of structures in both the ther- 
modynamic ensemble (a) and the mutational neighborhood 
(c) of the mdo-mir-1 miRNA precursor sequence have a sig- 
nificantly smaller fraction of structures that are highly dis- 
similar then sample sequences with identical MFE structure. 
Comparing the the averaged distribution of stem-loop (black 
bars) and random corresponding random sequence sample dis- 
tributions (white bars) in the mutational neighborhood (d) 
to similar averaged distributions in the thermodynamic en- 
sembles of the same sequences (b) shows that the tendency 
among miRNA precursor sequences for increased robustness 
is present both in the mutational neighborhood and the ther- 
modynamic ensemble, i.e. miRNA precursor sequences show 
excess robustness in the face of both thermal and mutational 
perturbation. 
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TABLE I: Phylogenetic breakdown of different measures of robustness 



group / species 


Tn Rn 


S n 


fa Rs 


S s 


ft (25) 7i t (25) S t (25) 


c(r s , n(25)) # of seqs. 


all 


0.44 0.59 


0.08 


0.29 0.78 


0.17 


0.31 


0.74 


0.28 


0.73 


3641 


vertebrate 


0.46 0.55 


0.06 


0.31 0.78 


0.13 


0.29 


0.75 


0.30 


0.76 


2215 


invertebrate 


0.37 0.68 


0.10 


0.21 0.88 


0.27 


0.22 


0.84 


0.36 


0.73 


488 


landplant 


0.41 0.63 


0.11 


0.31 0.75 


0.21 


0.40 


0.63 


0.19 


0.68 


848 


virus 


0.38 0.66 


0.09 


0.23 0.84 


0.18 


0.21 


0.85 


0.32 


0.65 


82 


Homo sapiens 


0.48 0.53 


0.04 


0.32 0.75 


0.11 


0.28 


0.76 


0.33 


0.74 


471 


Mus musculus 


0.46 0.56 


0.06 


0.33 0.76 


0.12 


0.31 


0.74 


0.27 


0.79 


373 


Drosophila melanogaster 


0.40 0.64 0.08 


0.22 0.88 


0.24 


0.23 


0.82 


0.35 


0.74 


78 


Caenorhabditis elegans 


0.30 0.78 


0.18 


0.20 0.89 


0.37 


0.23 


0.82 


0.34 


0.75 


114 


Arabidopsis thaliana 


0.39 0.62 


0.11 


0.29 0.78 


0.19 


0.43 


0.60 


0.15 


0.75 


131 


Epstein-Barr virus 


0.31 0.78 


0.00 


0.16 0.96 


0.22 


0.16 


0.87 


0.48 


0.81 


23 



Average rank-scores that indicate significantly increased according to both measures discussed in the Materials and Methods 

section (p- value < 10~ 3 ) are given in bold. 



sure the mutational robustness of a given sequence we 
use d the measures introduced by Borenstein and Rup- 
pin [Borenstein and Ruppin 2006j | : (i) the structural dis- 
tance based mutational robustness measure tj s of an RNA 
sequence of length L is defined by t] s = 1/ (3L) Y^to(L ~ 
di)/L, where d; is the base-pair distance between the 
secondary structure of mutant i and the native sequence 
(given by the number of base pairs present in one struc- 
ture, but not the other) , and the sum goes over all 3L sin- 
gle mutant neighbours and (ii) the more stringent mea- 
sure rj n is simply defined as the fraction of neutral single 
mutant neighbours, i.e. those that have identical MFE 
structure to the original sequence. In order to quantify 
the level of excess mutational robustness among miRNA 
precursor sequences we counted, for each miRNA pre- 
cursor sequence, the number of sample sequences that 
have higher mutational robustness according to a given 
measure (see Materials and Methods) and used this to 
calculate the rank scores r s and r„, defined as the frac- 
tion of sample sequences with identical or higher robust- 
ness according to rj s and r] n , respectively. To facilitate 
an overview of the extent of excess mutational robust- 
ness we also calculated the average of the rank scores 
over all miRNA precursor sequences f s and f n as well as 
the fraction of miRNA precursor sequences with higher 
than average robustness (i.e. rank-scores < 0.5) R s and 
R n and the fraction of sequences with statistically sig- 
nificant increased robustness (i.e. rank-scores < 0.05, see 
Materials and Methods) S s and S n , respectively, accord- 
ing to a give measure. The statistical significance of both 
rank scores for individual miRNA precursor sequences as 
well as that of the finding a given fraction of robust se- 
quences among a group of sequences was determined as 
detailed in the Materials and Methods section. 

Reexamining the mutational robustness of miRNA pre- 
cursor sequences in comparison to an unbiased sample of 



sequences with identical MFE structure we found that 
miRNA precursor sequences - in contrast to the results of 
Borenstein and Ruppin - do not have significantly more 
neutral single mutant neighbours than sample sequences, 
but do show a statistically significant increase in robust- 
ness measured according to T] s (see Fig.QJ) and Table UJ. 
In other words, native miRNA precursor sequences have 
on average the same number of single mutant neighbours 
with MFE structures identical to their own, as random 
sample sequences with the same structure. The MFE 
structure of those single mutant neighbours that are not 
identical to their own are, on the other hand, significantly 
more similar than in the case of sample sequences. 

The presence of excess mutational robustness is, by 
itself, insufficient to determine whether mutational ro- 
bustness has evolved as a result of direct selection or in 
congruence with selection for environmental robustness . 
As established previously Borenstein and Ruppin 2006j 
there is evidence for excess thermodynamic robustness, 
robustness to thermal fluctuations, as evidenced by a 
significantly lower than chance minimum folding energy 
among miRNA precursor sequences. Defining the envi- 
ronmental robustness measure t)e simply as minus the 
minimum folding energy we also find te = 0.278, Re = 
0.796, Se = 0.220 using unbiased sampling. The cor- 
relation between r s and te across miRNA precursor se- 
quences is, however, rather low with a Pearson's correla- 
tion coefficient of c(r s , te) — 0.217 and c(r n , r_g) = 0.071. 
The minimum folding energy is a somewhat crude mea- 
sure of thermodynamic robustness and does not even re- 
flect the excess mutational robustness according to the 
measure rj s . There is no good reason to assume that 
a low MFE in itself confers environmental robustness, as 
even for relatively high free energies a given sequence may 
none the less with high probability fold into structures 
sufficiently similar to the MFE structure to remain func- 
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FIG. 3: To quantify the extent of thermodynamic robustness 
present in miRNA precursor sequences we examined the rank 
statistics of the threshold thermodynamic robustness rjt(dth.) 
which we defined as the cumulative frequency of structures in 
the thermodynamic ensemble that are equal to or less than 
a distance threshold dth.. a For threshold values lager than 
d t h. ~ 20 the average rank of miRNA precursor sequences 
with respect to T]t(dth.)> denoted by ft (dth.), becomes larger 
than f s , while remaining highly correlated with it. For thresh- 
olds dth. > 25 the r?t(dth.) values for an increasingly larger 
fraction of random sample sequences becomes indistinguish- 
able from the r/t(dth.) value of the original miRNA precursor 
sequence due to a lack of structures with d > dth. in our finite 
resolution sample of the equilibrium ensemble. The fraction 
of sequences which remain distinguishable are labeled frac- 
tion resolved, b Among the majority of the random sample 
sequences that are resolved, however, the threshold robust- 
ness of large number of the miRNA precursor sequences be- 
comes highly significant (black and white bars) compared to 
the mutation robustness as measured by r/ s (grey bars) , whilst 
retaining a high correlation among r s and rt(d t h.). 



tional. The large number of miRNA precursor sequences 
that exhibit excess mutational robustness as measured 
by the structural similarity based measure rj s suggests 
that a strict adherence to the MFE structure is not nec- 
essary to retain functionality - a sufficiently similar, but 
not necessarily identical, secondary structure is enough 
to guarantee the excision of the proper subsequence. This 
is further supported by the fact that folding free energy 
alone is not sufficient to discriminate miRNA precursors, 
as well as recent evidence that a diverse set of structural 
features are needed for successful cleav age of a miRNA 
precursor sequence [Ritchie et al. 2007jj . 



To construct an appropriate measure of thermody- 
namic robustness that also reflects this observation we 
would need to know the extent of similarity that is re- 
quired to retain functionality - indeed we would require 
detailed knowledge of the interaction between the RNA 
substrate and the enzyme Dicer to establish an appropri- 
ate measure of structure similarity. As such information 
is not at present available we chose to use the most simple 
and widely employed structure similarity measure, the 
base-pair distance used above. In order to determine the 
extent of similarity required to retain functionality we de- 
fined the threshold thermodynamic robustness measure 
??t(dth.) as a function of the threshold distance dth., by 
equating it with the probability in the equilibrium ther- 
modynamic ensemble of structures that have base-pair 
distances equal to or less than a threshold dth. with re- 
spect to the MFE structure, i.e. 



(dth.) = $Z#(dth.-di) 



c -E % /kT 



(1) 



where the sum goes over the set of all possible structures 
il, di denotes the base-pair distance of structure i to the 
MFE structure, Z = X^esi c~ Ei ^ kT is the partition sum 
and H (x) is the unit step function, i.e. H (x) = if x < 
and H{x) = 1 if x > 0. 

Examining the thermodynamic robustness of miRNA 
precursor sequences in comparison to an unbiased sam- 
ple of sequences with identical MFE structure we found 
that miRNA precursor sequences have significantly more 
structures in their equilibrium thermodynamic ensemble 
that are similar to the MFE structure than sample se- 
quences (see Fig. [2k,b and Table UJ. In other words 
miRNA precursor sequences tend to adapt more similar 
structures as a result of thermal fluctuations than ran- 
dom sample sequences with the same structure. Calcu- 
lating the average rank score ft (dth.) and the fraction of 
robust Rt(dth.) and significantly robust St {dth.) miRNA 
precursor sequences, with respect to the measure r; t (dth.) 
(Fig. [3k) and examining the distribution of structures as 
a function of the base-pair distance for individual miRNA 
precursor sequences (see e.g. Fig.[2k) indicates that above 
a threshold distance d t h. ~ 20 the measures start to sat- 
urate, yielding an estimate of the required similarity to 
retain function. The correlation between the rank-score 
of miRNA precursor sequences according to the distance 
similarity based mutational robustness measure and the 
threshold thermodynamics measure is high for all thresh- 
old values. This is the direct result of the high degree 
of similarity between the distribution of structures in 
the thermodynamic ensemble and the mutational neigh- 
borhood (Fig. [2]). The average rank score ft (dth.) and 
the fraction of robust Rt(dth.) and significantly robust 
S t {d t h.) miRNA precursor sequences with respect to the 
threshold thermodynamic robustness measure indicate a 
markedly larger extent of excess robustness than their 
counterparts for mutation robustness, i.e. f s , R s and S s 
(see Fig. |3k,b and Table fl} above d t h. > 20. 
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Discussion 

The results presented above demonstrate the corre- 
lated presence of excess environmental (thermodynamic) 
and genetic (mutational) robustness among miRNA pre- 
cursor sequences as measured according to, respectively, 
r] s and r/t(d t h.)- A rather general causality between envi- 
ronmental and genetic robustness in the context of RNA 
secondar y structure has been sug gested by Ancel and 
Fontana [Ancel and Fontana 2000| . who studied the dy- 
namics of an in silico population of RNA sequences evolv- 
ing towards a predefined target shape. They found that a 
correlation exists between the set of shapes in the plastic 
repertoire of a sequence and the set of dominant (min- 
imum free energy) shapes in its genetic neighborhood. 
They argue that this statistical property of the RNA 
genotype-phenotype map, which they call plastogenetic 
congruence, traps populations in regions where most ge- 
netic variation is phenotypically neutral. In other words 
RNA sequences explore a similar repertoire of subopti- 
mal structures as a result of perturbations due to muta- 
tions and perturbations resulting from thermal fluctua- 
tions, and selection for a given target structure favours se- 
quences with higher robustness to perturbations of both 
type. 

Since, in contrast to genetic robustness, environmen- 
tal robustness does not require high values of uN e , as 
it is a property of the sequence and not its mutational 
neighborhood, we contend that the observed bias in mu- 
tational robustness is in fact the result of the congruent 



evolution of environmental and genetic robustness. 

The correlation between the response to heritable 
(mutational) and nonhcritablc (thermodynamic) pertur- 
bation, and hence the congruent evolution of genetic 
and environmental robustness may extend to other sys- 
tems with genotype-phenotype maps different from RNA 
secondary structure. In particular, Xia and Levitt 



Xia and Levitt 2002l | , have found compelling evidence of 



the correlated evolution of increased thermodynamic sta- 
bility and the number of neutral neighbours in lattice pro- 
tein models. Understanding the relationship between se- 
quence, structure, and function is, and will remain to be 
in the foreseeable future, a central theme in both molec- 
ular and evolutionary biology. A comprehensive view of 
how the relationship between sequence, structure, and 
function is shaped during the course of evolution must 
take into consideration both the potential correlations 
that arise from the physics of the structure-sequence re- 
lationship as well as the relevant population genetic con- 
ditions in the context of which it takes on the role of 
a genotype-phenotype map. In the context of compu- 
tational miRNA gene discovery our thermodynamic ro- 
bustness measure potentially offers an improved struc- 
tural feature that may perform better than the free en- 
ergy score of the hairpin or its ensemble diversity (which 
have proved uninformative Frevhult et al. 20051 ]). 
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